% function computes a vesselness filter on a widefield image of the brain
% and creates a rough mask to exclude the vessels

function [parenchymaMask,roughvesselmask] = getParenchymaPixels(im,windowmask,scaling,windowerosion)

if nargin < 3 || isempty(scaling)
    scaling = 5:45;
end
if nargin < 4 || isempty(windowerosion)
    windowerosion = 25;
end
maskedim = im;
maskedim(~windowmask) = prctile(im(:),80);
filtim = fibermetric(maskedim,scaling,'StructureSensitivity',25,'ObjectPolarity','dark');
roughvesselmask = imbinarize(filtim,'adaptive');

%% erode window mask to get rid of artifacts at edge:
SE = strel("disk",windowerosion);
windowmask = imerode(windowmask,SE);
roughvesselmask = roughvesselmask & windowmask;

%%dilate vessel mask just to be safe:
SE2 = strel("disk",5);

parenchymaMask = windowmask & ~imdilate(roughvesselmask,SE2);
